THE INVESTIGATION OF OPTimL DISCRETE APPROX I »MT IONS 
FOR REAL TIME FLIGHT SIMULATIONS 

Final Technical F^>ort 
Grant No. NASA ^SG 1151 

Submitted to; 

NASA Scientific 4 Technical Information Facility 
P. 0. Box 8757 

Baltimore/Washington International Airport 
Maryland 21240 

Submitted by: 

E. A. Parrish 
E. S. McVey 
G. Cocrfc 

I K. Henderson 

T81 IIVISTIGITIO* OF *76-19158 

OPTIliL DISCBMl APPtCIIBATlOKS FOB BB&L 
TIBB FLIGHT SIBDLATIOIS Flcal Technical 

Beport (Ticginia onie.) 80 p BC $5*00 Onclas 

CSCL OIC G3/08 20654 

SCHOOL OF ENGINEERING AND 
APPLIED SCIENCE 


RESEARCH LABORATORIES FOR THE ENGINEERING SCIENCES 


UNIVERSITY OF VIRGINIA 


CHARLOTTESVILLE, VIRGINIA 22901 


Report No. EE-4041 -1 02-76 
March 1976 


THE INVESTIGATION OF OPTIMAL DISCRETE APPROXIMATIONS 
FOR REAL TIME FLIGHT SIMULATIONS 

Final Technical Report 
Grant No. NASA NSG 1151 

Submitted to: 

NASA Scientific & Technical Information Facility 
P. 0. Box 8757 

Baltimore/Washington International Airport 
Maryland 21240 

Submitted by: 

E. A. Parrish 
E. S. McVey 
G. Cook 
K. Henderson 


Department of Electrical Engineering 
RESEARCH LABORATORIES FOR THE ENGINEERING SCIENCES 
SCHOOL OF ENGINEERING AND APPLIED SCIENCE 
UNIVERSITY OF VIRGINIA 
CHARLOTTESVILLE, VIRGINIA 


Report No. EE-404 1 -1 02-76 
March 1976 


Copy No 


/ 



I. INTRODUCTION 


This report summarizes the results obtained during the first year 
of a continuing effort on the general topic of discrete approximations 
for real-time flight simulation. The report is divided into five Riajor 
topics as follows: 

1. Digital Autopilot Modelling — consideration of the particular 
problem of approximation of continuous autopilots by digital autopilots. 

2. Frequency Domain Synthesis of Discrete Representations — use of 
Bode plots and synthesis of transfer functions by asymptotic fits in 

a warped frequency domain. 

3. Substitutional Methods — an investigation of the various sub- 
stitution formulas, including the effects of nonlinearities. 

4. Use of Fade approximation to the solution of the matrix expo- 
nential arising from the discrete state equations. 

5. An Analytical Integration of the State Equation Using Interpo- 
lated Input — uses polynomial approximations to the input signal and 
integrates the state equations analytically. 

Each of the sections is self-contained in that the developments 
and conclusions are contained in each section. 

The work presented In this report represents the initial efforts 
and preliminary Investigations of the topics covered. The actual appli- 
cation of the various techniques has been limited to some rather simple 
linear and nonlinear systems, while the emphasis has been placed on 
theoretical investigations. This emphasis will shift somewhat in future 
efforts as more complex and nonlinear systems are simulated with the 
techniques which appear most promising as a result of the work presented 
in this report. 


2 



II. DIGITAL AUTOPILOT DESI»I 


Introduction 

The purposD of this part of the research Is to establish criteria 
for approximating continuous autopilots with digital autopilots. ,t 
is desired that the performance of an aircraft with a digital autopilot 
be similar erx>ugh to the continuous autopilot which it is to replace 
that pilots cannot perceive a difference between them. This is the 
first of three levels of work being considered in the investigation of 
digital autopilots. The three levels are: 

1. The problem of replacing existing continuous autopilots with 
a digitally implemented autopilot that is as nearly indistinguishable 
as practical from the continuous autopilot. 

2. Design of autopilots using digital control methods based on 
original performance specifications. 

5, 5rA+ir>n coptro! s'/stems and use of modern 

methods to control and optimize systan performance for flight situations, 
Sjch as optimal fuel control which existing autopilots do not provide. 

Although the ideal situation is for the performance of new digital 
autopilots to be identical to the continuous systems they replace, it is 
easy to prove that, in general, it is not possible to find exact dis- 
crete equivalents to continuous systems. The principal differences in 
the output signals will be: 

1. phase shift and attenuation due to the zero order hold at the 
computer output. 

2. phase shift and gain due to the digital representation of the 
continuous autopilot transfer function. 

3. aliasing of high frequency input components. 

4. distortion components due to sample-hold operations. 

Although these can result in major differences in the results of 
the two implementations. It should be possible to tnake them conform to 
commonly used control specifications to within any accuracy desired if 
the sampling frequency of the digital system is made high enough; but. 
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unfortunately, this parameter is usually restricted because it is desired 
to keep the scmpling frequency as low as practical to minimize computer 
capacity. Sample frequency Is a major parameter in the study. 

A logical approach and the one used here is to use design specifi- 
cations of continuous autopilot systems, but keep In mind also that the 
origiral design was based on continuous control technology and should 
take into account any new factors Introduced by a digital controller. 

The original specifications are not available but they undoubtedly C2.l!] 
consisted of such factors as phase rargin, gain margin, magnification, 
closed loop pole locations, rise time, overshoot, delay time, settling 
time, mean squared error to stochastic inputs, final value of error, 
dynamic tracking error, and bandwidth. 

investigation of Gain and Phase Shift Errors Due to Discretization 

Discrete approximation techniques to represent continuous functions 
in digital computers are studied and compared in the following. As 
already stated, the objective is to design digital controllers which will 
substitute for existing analog controllers, such that the digital con- 
trollers are as similar in function as practical to the hardware they 
replace. Evaluation on an input-output basis will make use of phase 
and gain specifications. 

Consider the autopilot transfer function for a commercial aircraft 
as shown in Fig. 2.1. It is 



36(s + l.65)(s^ + 2.31s + 2.72) 

(s + 0.62)(s2 + 5.62s + 3.l)(s2 + 8.4s + 36) 


( 2 . 1 ) 


where U(s) is the Laplace transform of the pitch rate input signal and 
Y(s) is the Laplace transform of the control signal. The Bode diagram 
for the above transfer function is shown in Fig. 2.2. It is inferred 
from the figure that the crossover frequency for the closed loop system 
should be located at about id = 3 rad/sec. The crossover frequency is 
of special importance and will be used in later discussions. 
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Ftgure 2.1 Autopilot Transfer Function for a Commercial Aircraft 
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Discrete Approximation 


If G(s) is a continuous transfer function as shcMfn in Fig. 2 
its discrete counterpart would appear as in Fig. 2.3(b) where 

Y(s) = G (s)E(s) 
c . 

Y»(s) = CG^(s)E(s)D« 

Y.»(s) = G «(s)E»(s) 
d A 

The ideal discrete system for fhe applications here would be 
that yields 

y'(t) = y(t) for all t > 0 

However due to the discrete nature of the digital controller 
It is impossible to realize Eq. (2.5) exactly. The best possible 
Crete system could give, instead, 

y’(t) = y(t) for t = 0,T,2T, ... 

where T is the sampling period. Because 

y'(nT) = y^j(nT) n = 0,1,2, ... 

in this case Eq. (2.6) is satisfied 

y(nT) = y^jtnT) n = 0,1,2, ... 

or 

Y»(s) = Y^»(s) 

Using Eqs. (2.3), (2.4), and (2.9), 

CG^(s)E(s) 3« = G^«(s)E*(s). 

Taking the z-transform of both sides of Eq. (2.10), 

Z{G^(s)E(s)} = G.(z)E(z) 

C n 


3(a), 

( 2 . 2 ) 

(2.3) 

(2.4) 
one 

(2.5) 
dis- 

( 2 . 6 ) 

(2.7) 

( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 
( 2 . 11 ) 
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Figure 2.3 (a) Continuous Transfer Function and (b) its Discrete Counterpart 




or 

Z{6 (s)E(s)} 

\(7) ^2.12) 

where Z{*} indicates the z-transform. 

In Eq. (2.12), D(z) cannot be obtained independent of the input 
function e(t). This justifies the previous assertion that an exact 
digital equivalent of the continuous autopilot cannot be obtained. Thus 
one must settle for a difference equation transfer function that gives 
an acceptable approximation to the solutions for an arbitrary input, 
eft). 

Various discrete approximation techniques have been proposed to 
transform a previously designed continuous controller into a discrete 
equivalent. These discrete approximation techniques can be divided 
into three broad categories: (I) numerical methods (2) operational 
methods (3) input approximation methods. 

Numerical methods generally provide a means of obtainihg accurate 
approximation. However, since these methods generally take a great deal 
of calculation time, their application to discrete approximations is 
often limited. In the operational methods, every integrating operator 
of the continuous transfer function is replaced by a discrete integrating 
operator in order to obtain the discrete transfer function. In the 
input approximation methods, the input is assumed to be approximated 
in a certain way, for example, by a stair-step function, or by straight 
line segments, etc. 

In the literature, Tustin’s method is generally preferred among 
the operational methods and the linear segment input approximation method 
is considered to be one of the best input approximation methods. These 
two methods will be used in the following discussion to show that the 
degree to which the digital control system approximates the continuous 
system is determined primarily by the zero-order hold rather than by 
the discrete representation (i.e., discretization) of the continuous 
controller transfer function. Various methods for discretizations are 
compared later. 
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Tustin's Approximation 

Tusttn's method enjoys the merit of simplicity. This method is 
usually quite accurate, does not Introduce spurious solutions, and is 
ideally suited for operational calculus operations due to Its cascading 
property. To apply Tustin's method, first divide the numerator and 
denominator by the highest power of s in the denominator, s^, then 
replace < ^ by ( )* to obtain the discrete transfer function. 


For G(s) consisting of a simple pole at -a. 



Straight Line Approximation 

This technique assumes that the excitation can be approximated by 
a series of straight-line segments as shown in Fig. 2.4(a) for an 
arbitrary function. The piecewise linear input may be applied to any 
system by placing a sample and a triangular hold (unrealizable first 
order hold) before the normal input as shown in Fig. 2.4(b). For the 
above example 
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_ (r - D- 

TT~ 


(z » I , Ca(z - e~^^)Tz - z(z - IHz - e + z(z - 1)^] 
Tz a2(z - l)2(z - 


aT(z - e~^^) - (z - l)(z - e"®^) + (z - 1)^ 


Ta^(z - e~®') 


(e"®^ + aT - l)z + (I - e"®^ - aTe"®^) 


Ta^(z - e 


(2.14) 


for the stratght-l Ine approximation. 

Gain and Phase Shift 

An approximation of continuous autopilots with digital autopilots 
should Include matching the gain, phase shift, and steady state errors 
of the systems. 

The gain and phase shift of the digital and analog controllers In 
the frequency domain will now be considered. 

Referring to Fig. 2.3(b) and using Eq. (2.4), 


Y'(s, = G. (s)Y.*(s) 
ho d 


(2.15) 


= G^^(s)G^*(s)E*(s) 


(2.16) 


Y'(j(o) = Gj^^(Ja))G^»(Ja>)E*(Jaj) 


The frequency range of interest Is 


(2.17) 


U < Cl)g/2 

where u Is the input frequency and Ug Is the sampling frequency. 
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If the sanpling theorem Is satisfied (and It will be for cases of 
Interest here) 

E»(j«) £(j«) for « < e^/2 (2.18) 

Thus Eq. (2.17) becomes 

Y*(j«) = 6j^(j«)6^*(j») YE(jm) (2.19) 

or 

Q 

1“ (j«) = (je)G^*(jw) (2.20) 

For the continuous transfer function in Fig. 2.3(a) 

Y 

^ (je) = G(j«) (2.21) 

The transfer function obtained from Eq. (2.20) is to be matched as 
closely as possible to the transfer function of Eq. (2.21). 

For the gain 

Y» l^iw.^j“^l 

I U«)| = . ’f |G^»(jM)| (2.22) 

and 

1 p = |G(j«)j (2.23) 

Therefore 


ly^i r— -i-BTisr 

For the phase shift 


'« < «g/2 


I (Jm) = /G^o(J<a)/T 



(2.24) 


(2.25) 
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and 


I ^ (jg» = I G(ju) 

Therefore 

(>' - I g <>) 

w < Wg/2 

As indicated in Eq. (2.24) and Eq. (2.27) the differences in the gain 
and the phase shift of the digital controller and the analog controller 
can be divided into two terms; one is due to G. (ju); the other is due 
to G^*( ju)/G( ju) . The former depends only on the sampling frequency, 
whereas the latter would depend on the specific continuous transfer 
function to be discretized and on the specific discretization technique 
used as well as on the sampling frequency. It is interesting to compare 
the differences due to the two terms above in a specific example. 

Consider the transfer function 


(2.26) 


+ j (G^*(j(o) - G(jg))) , 
(2.27) 


G(s) 


I 

s + I 


Tustin's approximation from Eq. (2.13) with a = i is 


Gat(z) 


T(z + I) 

T(z + I) + 2(z - I) 


For fg = 5, T = .2 

,, _ 0.2(z +1) _ 0.090909(z + I) 

~ 0.2(z + I) + 2(z - \) z - O.SlSl82 


G^*( ju) 


0^090909_(eJ°'^“ + |) 
JC.2oj 


0.818182 


(2.28) 


14 



For fg « 10. T = .1 


0.047619(0 


®AT ‘J“’ — j37nr 


JO.I. ^ 


I) 


- 0.904762 


(2.29) 


For * 20, T * 0.05 


6AT*(j ) = 


0.02459(6“^°*®^ + I) 

J6.b$ 


- 0.95122 


(2.30) 


The straight-line approximation from Eq. (2.14) with a = I is 




(T - I e~^)2 t I - (T t I )e~^ 
T(z - e"^) 


For * 5, T = .2 

9 


^ _ 0.0187312 + 0.017523 _ 0.095655z + 0.087615 

“ 0.2(z - 0.818711 > z“-ir.675751 


GAS*(j«) = 


0.0936p(e-^°*^“ + 0,935508 
eJ0-2“ _ 0.818731 


(2.31) 


For = 10, T = .1 


GAs*(j«) 


0.048^(6-^°^'“ + 0.96734) 
. 0.904837 


For fg = 20, T = 0.05 


(2.32) 


GAs*(ju.) 


0.02458(6-*°’°^“ + 0.98373) 
eJ0*05u) _ 0.951229 


(2.35) 


For the zero-order hold 



< ju») 


2n .-in(«/..) 

b) IKw/w ) 
s s 
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(2.34) 


sin n(u/(i) ) 

|6ho<J»>l = T • I lfU.^1 I 

and 

/6^(jw ) = -n( ^ > sgnCsin H{a/Ug>D (2.35) 

And, from Eq. (2.34) 

“T 

For the frequencies of interest, 
sgnCsin IKw/m^)] = I. 

Therefore, 


I sin n(a/a„) 

I _ s 

n(»/Wg) 


(2.36) 


L 


G. (ju) = -n( -52- ) 
ho ** Wg 


(2.37) 


The results are shown in Figs. 2.5, 2.6, 2.7, 2.8, and 2.9. Fig. 2.5 
shows |G(Ju) I and lG^*(ju)| as a function of frequency for different 
sampling frequencies. The distinction between the two different approxi- 
mation methods has been ignored because the differences are negligible for 
the present purpose. This is discussed later in more detail. It is noted 
that in the low frequency range, |G^*(jw)| follows |6(ju)| almost perfectly 
and in the high frequency range, |G.*(ju)| fol lows |G(ju)| more closely as 
the sampling rates get higher as expected. Figure 2.6 shows | ho| and 

I 

jG^*(j(i))/G(jw) I as a function of frequency. From Fig. 2.7 to Fig. 2.9 
yG^^(jM) and ^G^*( Jm) - ^G( Jm) are shown for different sampling 

frequencies. It is noted that j G^^(ju) is much larger than 

/G/(jg) ) - ^(Jm) over most of the frequency range for all sampling 

frequencies. 
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Figure 2.6 Attenuation Due to Discretization Compared with that due to Zero-Cr:er r-cld. 

Tustin and Straight-Line ADcroxfratFon Methods are used for the Discretization 





Figure 2,7 Phase Shift due to Discretization Conpared with +hat due to Zero-Crder t'cl: ■‘or 
fg = 5 samples/sec 




Figure 2.S Phase Shift due to Discretization Compared with that due to Zero-Order ifold 
for * 1C samples/sec 



RKPROl)i't--‘- OF THE 
fS POOR 



C.OOI 0,01 0.1 I 10 too 


♦ u (rad/sec.) 

Figure 2.9 Phase Shift due to Discretization Compared with that due to Zero-Order Hoid 
for fg = 20 sampies/sec 



Consideration of Figs. (2. 5-2. 9) shows that In this discretization 
process the phase shift would be a more important determining factor 

I 

than the gain. For example, for = 5 samples/sec, at w = 3 rad/sec 

/e,^(j ) . -17.3° 


and 


y G^*(jtti) - / G(Jtt)) = 0® (stralght-l Ine app.) 

= -0.5® (Tustin's app.) 



0.99 


|G^*(ja»)/G(ju)| * 0.975 


Thus 


and 


/ShoU"> * / V<J“> 


/g(w) 


-17.3® 


-^r — 


|G/(jo)) 

I G'Oo)l 


0.97 


The change in phase shift by -17.3® would produce more significant effects 
than the change in gain by 0.97 in most practical systems. In fact, it 
is found that over most of the frequency range of interest the phase 
shift Is more important than the gain. 

It is important to note that in the range of useful frequencies the 
change in phase shift due to discretization comes from the zero order 
hold and one cannot choose between the above two approximation techniques 
on the basis of performance. A designer would be justified in using the 
simpler of the techniques. 

It can also be demonstrated for general cases using Tustin’s approxi- 
mation that the value of ^G^*(j(«)) - /G(jto) is negligible compared to 

/ ^ho^*^“^ * demonstrarlon proceeds as follows. 
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Assume a transfer function with a simple pole. 


G(s) = 


s + a 


/ G(.ih») = - tan-l J 


(2.38) 


Then, using Tustin's method 


Ga(z) = 


T z + 1 

7 


7TT T z+T 

' ^ ® * I ‘ 


T(z + I) 

5(z - I i + at(z + 1 ) 


T(eJ“^ + I ) 


G *(jo)) = G.(z) . , = — ’-T 

^ '' z » eJ“^ 2 (eJ“^ - I) + aT(eJ“^ + |) 


e' 


JuVi -ja>T/2 


- e 


;W77—=W72 


+ aT 


2jtan( ) + aT 


(2.39) 


The phase angle of the approximation is given by 

G/(jio) = -tan-^C ^ tan ( ^^ )] (2.40) 

Now, Eq. (2.40) should be compared wtih Eq. (2.38). Let 0, and 02 be 
defined as 

0j = ~tan“^ ^ ©2 “ “tan ^ E tjy tan( )] (2.41 ) 


Then 


tan ©2 ” tanj 

tan(02 ©i) I + tan ©2+an ©i 

_ - ^ tan ( !^ ) . ( I ) 

' * < ?r f ’< ? ’ 


(2.42) 


23 


P.KPRODIICjr^'^' OF TBF 

'V- ■ ■ '■ 



Now show that the numerator of Eq. (2.42) Is negligible. First, 


s 

= n( -2L Hrad) 

“s 

= -J- X 180® (2.45) 

Referring to Eq. (2.57), wT/2 is the phase lag due to the zero-order 

hold. 

If ( ^ ) --5® (which should be done to realize a digital controller 
that approximates the continuous case), the tangent of the angle may be 
replaced by the angle so 


tan ~ = (oT/2 

and from the numerator of Eq. (2.42) 


- 2 if ) + £> 

aT 2 a 


2 u)T . (D 


a a ~ 


0 


I .e., 01 - 02 = 0 
or 

01 = 02 


(2.44) 


For the case where the transfer function has a pair of complex 
poles, a similar development is possible. 

Because Tustin's approximation has a cascading property, i ,e. 

^Tust:,i^^l^^^^2(s)] = {Zyy^^j^CGi(s)3HZ^^^^j^[G2(s)]} ( 

the above approximation holds for general transfer functions. 
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Further Investigation of gain and phase differences Introduced by 
the discretization process for the example of the simple first order 
pole Is summarized in Figs. 2.10 through 2.16. 

Figure 2.10 shows the ratio jo))/G(Jii>) | for several values of 
system time constant as a function of the ratio of frequency to sampling 
frequency. Tustin's method was used for the approximation. The normalized 
magnitude of the zero-order hold is also shown for comparison. 

Figure 2.11 is a plot of the phase for the same conditions. This 
figure reiterates the conclusion that the phase lag Introduced by the 
zero-order hold predominates. 

Figures 2.12 and 2.13 differ from Figs. 2.10 and 2.11 only In that 
the linear segment input approximation was used for the system. Again 
the phase lag of the zero-order hold predominates. 

Figures 2.14 through 2.15 compare the gain and phase shift of 
several discretization methods for two sampling frequencies. These 
figures show that one of the simpler methods (Tustin) is also one of 
the more accurate methods considered for this application. 

Selection of the Simulation Increment (Computer Sample Period) 

One of the most Important criteria in the choice of a particula.- 
digital simulation method is the length of the simulation Increment re- 
quired to produce a given accuracy. It has been well documented in the 
literature and in our efforts that some methods require a smaller simulation 
increment than others to produce equivalent results, and , in fact, some 
methods are unstable for increments that produce good results in other 
methods. Therefore the choice of the smallest simulation increment to 
produce a given accuracy has been the subject of recent Investigation. 

This topic is of particular importance when the simulation must be 
performed in real time and on the smallest possible machine. 

The structure of the problem of finding the smallest possible simula- 
tion increment which will yield a specified accuracy for a given simulation 
method suggests the possibility of using classic optimization theory. This 
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Figure 2.11 Phase Difference for the Same Conditions as *^lgure 2.i0 
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Figure 2.13 Phase Difference ♦or Sane Conditions as Figure 2.12 






BEPBQDUCIBlLnY OP THB 
0BIQ1NAL PAGE IS POOR 



Figure 2.14 p^ 1 ase difference due to Discretization for Several Methods (Sampling frecuenoy of 5 Hz) 
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Figure 2.16 Attenuation due to Dlscretlzatl 
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was the approach taken by Rosko and Darling C2.23 and Rosko C2.33. Their 
results are currently being impieR»nted to aid in the choice of simulation 
method and for possible inclusion in the Simulation Design Package. 

Rosko and Ourling's method is based upon the assun^tion that round*^ff 
errors have been or can be made negligible for the time increnent to be 
determined. The problem then is to find a compromise between o»> 4 >utation 
time and simulation error. The method will now be summarized. 

Formulate a performance index which will mathematically reflect the 
parameters to be minimzed and their relative in 4 X>rtance. Functionally, 

J = f(E,T, ...) (2.46) 

where E is a measure of the simulation error and T is the simulation in- 
crement to be determined. An example of such a perfomance measure is 

J = E + (2.47) 


where 3 is a constant weight. 

For deterministic inputs the integral squared error for a system 
simulated with interpolated intervals is given by C2.2, 2.3] 

E(t) = e^(t) = /^Cy.(t) - y^(t)]2dt (2.48) 

where Vj(t) and y^(t) are the outputs of the ideal and approximated 
systems respectively. Using Parseval’s theorem we may express the inte- 
gral squared error as 

E(t) = ^ iZo <2.49) 

where 4^(i«>) is the spectral density of the error. 

For the case of discrete simulation where the error Is considered 
only at the sampling intervals, the summation of squared error is given 
by [2.2, 2.3] 

* I 

E(nT) = J e2(nT) = -s— ^_E(z)E(z”< )z"*dz (2.50) 

n=0 J 
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which may bo evaluated by sunmlng the residues of the poles of the 
Integrand inside the unit circle in the z-plane. 

Select a reasonable increment value, T = and calculate the 
initial value of the error using the appropriate error criterion. This 
gives an initial value of the performance index. 

Jl = El + 1^ . (2.51) 

(^iculate an estimated value of the error and performance index by 

E 2 = y(T) = ^ Ell (2.52) 

* |T=T2 

J 2 = E 2 + Y I <2*53) 

|t=T2 

A necessary condition for optimality is 

'*1 a T B 1 


which yields 



(2.55) 


Use T = T 2 to predict E using the appropriate error formulation. 
Then calculate 


J 2 = E2+ ^ 


subject to the constraint 


AJ 2 

J 2 


J2 - Ji 


i To 


(2.56) 


(2.57) 


where Tq is a constant to be chosen. 

For all succeeding steps, m > 3, the following procedure applies: 
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E » V(T) = 
ni 


}S, •" ‘ 


(2.58) 


where Lagranglan Interpolation has been used, and 


n = KT-TiHT-To) ••• (T-T ) 
m m 


(2.59) 




(2.60) 


2) Now calculate 


A A O 

Jn. = + T- 

m mi. 


(2.61) 


3) Calculate using Newton-Raphson quasi I Inearization 

where 

J'(T_ ,) 

T , = T - for k = I (2.62) 

m, I m-i j • • (T ) 

m-l 


J'(T , ,) 

T ^ = T for k > I 

m,k m,k— I j 

m,k-l 


(2.63) 


subject to the constraint 


T , - T , , 
m,k m.k-l 


(2.64) 


where tq Is a chosen constant. 

4) Calculate the predicted value E using T = T and the appropriate 

m m 

error formulation. 

5) Calculate the performance Index 


J = E + ^ 
m m T 

m 


(2.65) 
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subject to 




1 r© 


( 2 . 66 ) 


The calculation terminates whenever Eq. (2.66) Is satisfied. 


It has been shown that the discretization method is unimportant 
from a performance view based on phase shift and gain considerations. 

This conclusion will also be checked using the procedure presented above. 


Ccnclusions and Vtork in Progress 

The performance of a continuous autopilot which is implemented 
digitally is affected so much by the zero order hold that the discreti** 
zation method is a secondary consideral ion. Tustin is a relatively simple 
method that is satisfactory. Computer speed (sampling rate) may be 
established on the basis of the phase shift a designer will allow to be 
introduced by the zero order hold. Phase lag decreases as sample speed 
increases. Gain constants should be maintained to keep steady state 
errors constant. The relative stability of the aircraft will be de- 
creased by the digital implementation so the value of computer sampling 
time should be based on the decrease in phase margins a designer is 
willing to allow. This value can be determined through sensitivity 
studies. Work is beginning on an example sensitivity simulation to 
demonstrate how one can approach •'■he problem of specifying an allowable 
change in phase margins. Optimum simulation increment work is also 
underway to give credence 'lo the above conclusions from another point 
of view. Round off error work Is to be done during the next year. 

Further studies of discretization methods considered in this section 
will be undertaken for more complex (higher order) systems with a 
view toward possibly utilizing phase lead introduced by some of the 
methods to offset phase lag Introduced by the zero-order hold. 
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III. FREQUENCV DOMAIN SYNTHESIS OF DISCRETE REPRESENTATIONS 

Analysis of linear, continuous time systems Is frequently done In 
the frequency domain. By use of the Laplace Transform, differential 
equations are replaced by algebraic equations which. In general, are 
simpler to solve. In designing a discrete time system to approximate 
the performance of a particular continuous time system. It would be 
helpful If some of the familiar analysis techniques could be utilized. 

A design procedure Is presented by which a discrete time transfer 
function can be developed In the frequency domain. The resulting system 
will have frequency domain characteristics similar to the continuous 
time system from O.C. to one half the sampling frequency. It will also 
be shown that the time domain performance of the two systems will be 
similar. The chaiacteristlcs of the type of continuous time system 
which can be most closely approximated using this design method will 
also be discussed. 

A simple design example will be presented to explain the methodology 
of the design procedure. Following that, the frequency domain design 
and time domain evaluation for the autopilot will be discussed. 

Explanation of Design Procedure 

The transfer function for this example is shown In Figure 3.1 

F(s)=j4^ (3.1) 

and the sampling frequency Is chosen (arbitrarily) to be 20 rad/sec. 

The design objective Is to develop a discrete time transfer function 
F(z), whose frequency domain character' sties will closely approximate 
those of F(s)from0 to 10 rad/sec. A graphical version of the design 
procedure will be presented to Illustrate the method. Part or all of 
the procedure can be automated. 

The first step Is to plot the magnitude of the transfer function 
versus frequency, i.e., 20 log F(Jw) versus log oi. The plot need not 
extend higher in frequency than one half the sampling frequency and 
may go as low In frequency as desired. 
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The next step is to define a complex variable p such that 


P = jX * (3.2) 

From Z- transform theory 

z = exp(sT) 

and for real frequencies this becomes 

z = exp(jwT) (3.3) 

where T = Zir/ai^ Is the sampling period, and Is the sampling frequency 
In radians per second. Substituting Eq. (3.3) Into Eq. (3.2) yietds 

p = jX = (exp(jtoT) - 1 )/(exp(ja»T) + I) (3.4) 

3y factoring and utilizing trigonometric identities, Eq. (3.4) can be 
expressed as 


p = jX = jslr.(a>T/2)/cos((oT/2) = jtan(wT/2) (5.5) 

The magnitude plot of Eq, (5.1) Is now re-plotted versus the fre- 
quency X s’-Ch that the following relation Is true (see Figure 3.2); 

F(jX)| = F(jo)T) (3.6) 

1 X- \j=tanWjT/2 

In other words the magnitude of F(juj) at a particular frequency u. is 
mapped Into a point on the F(jX) plot, having the same magnitude and 
occurring at a frequency Xj, where 

Xj = tan(w.T/2) (3.7) 

When 0 ) = 0, X = tan'O) = 0, and when in = tn^/2, X = tan(n/2) = <». Thus, 
the frequency ra-je 0 < u < m^/2 maps into the frequency range 0 < X < ®. 

The next step is to synthesize an equation for F(p) by inspecting the 
magnitude plot just made. The plot of F(jX) versus X will look similar 
to the plot of F(joj) versus u at low frequencies. As X increases, F(jX) 
will approach the value that F(Ju)) has at m^/2 along a horizontrl asymptote. 
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Equation (3.9) can be written in the form 

z + K 2 

F(z) = Ki -- n r- (3.10) 

* z + K 3 

where = (Xi/X 2 )(X 2 + l)/(Xj + I), K 2 = (X 2 - l)/(X 2 + I), 
and K 3 = (X| - l)/(Xi + 1). 

The frequency response of F(z) can be determined by setting Z = expijojT) = 
cos(uT) + jsin(wT). This results in the following two expressions: 

nKjgCF(z)] = 20 log{KiC(cos«T + K 2 )^ + (sinuT)2]’S/|:(cosMT + K 3 )^ + (slnuiT)2]*S) 
angrF(z)H = trn"^Csin(oT/(cosuT + K 2 )D - tan'^CsinuT/icosuT + K 3 )] 


Since z = exp(juT), these iwo functions can be plotted versus the original 
frequency u. The magnitude plot of F(z) should be quite close to the 
magnitude plot of F(s) from DC to w^/Z. Since the numerator and denominator 
of F(z) are of the same order in z, the phase shift will be zero degrees 
at 0)^/2. To try and shape the phase shift of F(z) to more closely 
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match F(s)» an all>pass filter can be added to F(p). This results in 
the n^ transfer function 


F(p) 


Xl P + X2 X3 - p 

< — TTT- > < y—r- 
X2 P + Xl X3 + p 


(3.11) 


making the substitution p - (z - l)/(z I) into Eq. (3.11) gives the 
fol i<Mlng: 


(z + K2 )(z + K3) 
(z + K4 )(z + K5) 


(3.12) 


where 


Ki = (Xi/X2)(X2 + l)(X3 - |)/(Xi + l)(X3 + I), 

K2 = (X2 - l)/(X2 + I), 

K3 = (X3 + l)/(X3 - I), 

= (Xl - l)/(Xi + I), 

and 

K5 = (X3 - l)/(X3 + I ) = l/K . 

The phase of this new F(z) will go to -180 degrees at w = a^/ 2 . The 
actual phase characteristics of the original F(s) will determine how 
high In frajuency the phase of F(z) matches that of F(s). Figure 3.3 
shows the magnitude of F(z), as well as the phase with and without an 
a 1 1 -pass f I Iter. 

To evaluate the time domain performance, F(z) can be expressed as 
a ratio of polynomials 

F(z) = KiCz^ + (K2 + K3 )z+ K 2K3yCz2 + (K^ + Kjiz + K4K5D 

(3.13) 

Multiplying the numerator and denominator of Eq. (3.13) by z“^ yields 


F(z) = = Kid + (K2 + K3 )z“^ + K2K3Z~^/D + (K4 + + K4K5Z-2 

(3.14) 


where C(z) and R(z) represent the output and input of the transfer 
function, respectively. Cross multiplying Eq. (3.14) produces the following 
expression. 
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D + (K4 + K5)2~^ + K^Ksz'^lCCz) = KiCl + (K 2 + Kj)z-‘ + KjKjZ-^lRCz) 

Applying the Real Translation Theorem from Z-Transform theory yields 
the foll^ing time domain recursive equations: 


CCnT] + (K + K )CC(n - I )T3 + K K CC(n - 2)T3 

= KiCR(nT) + (K 2 K3>RC(n - I )T] + KaKsRCCn - 2)T] 

(3.16) 

2 2 

CCnT] = - I a.CC(n - i )T] + Kj I b.R[(n - i )T] (3.17) 

i=l ' i=0 ' 


where the a's and b's are obvious from inspection of Eq. (3.16). By 
means of Eq. (3.17), the output C(nT) can be found for any input by 
merely specifying the sequence R(nT) as the integer n varies over the 
range of interest. 

The location of the poles ,nd zeros of F(p) and the location of the 
all-pass filter break frequency can be optimized, using whatever per- 
formance criterion desired. 


Application to Autopilot 

The design procedure outlined on the previous pages will now be 
applied to the pitch portion of the continuous time autopilot. The 
open- loop transfer function for this system Is 

H(s) = 

36 s^ + 2.31s + 2.72 s + I .65 s" + 7.25s + 81 

s2 + 8.4s +36 s2 + 5.62s +3.1 s + .62 I.IZSs^ + 13.33s + 81 

(5.18) 

The sampling frequency was chosen to be 40 rad/sec. This allows the range 
of interest to extend to approximately twice the highest critical fre- 
quency of the autopilot. 
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The fictitious transfer function H(p) was synthesized in the 
fol lowing manner: 

1. magCH(ju)] for 0 < « < 20 was calculated, as well as the 
frequency X, using the relation X = tan(cijT/2). 

2 . A plot of magCH(jX)] versus X was made (Fig. 3.4). 

3. Real poles and zeros of H(s) at u. were transformed into real 
poles and zeros of H(p) at X., where X. is given by Eq. (3.7). 

4. Complex roots of H(s) with undamped natural frequencies o< 

were transformed into complex roots of H(p) with natural frequencies also 
given by Eq. (3.7). The values of the damping ratios were preserved in 
synthesizing H(p). 

5. An initial value for a second order real zero was made by 
Inspection of the magnitude plot made in step 2. This additional term 
is due to the asymptotic approach of H(p) to the value of H(s) at 

(0 = u^/2. 

6. The substitution p= (z- l)/(z+ I) was made in the expression 
for H(p). 

7. The magnitude and phase plots of H(z) versus w are made. 

8. An initial choice for the all-pass filter break frequency was 
made, and the phase curve re-plotted. 

At this stage the magnitude curve of H(z) differs from that of H(s) 
by a maximum of 1.7 dB at a frequency of 6.5 rad/sec. The phase curve of 
H(z) has the same general shape as that of H(s) up to approximately 15 
rad/sec, but has an error of about 17-22 degrees from 8 to 15 rad/sec. 

It was decided at this point to optimize several parameters of H(z). 

The locations of the second order real zero added during the design 
process and the term resulting from the single real pole of H(s) it 
ui = 5 rad/sec were simultaneously optimized. The performance criterion 
used was the mean absolute error between the frequency domain magnifude 
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of H(s) and the magnitude of H(z) over the frequency range 0 to u^/2. 

The values of these two parameters were independently stepped through 
the range 0 to 20 in intervals of O.^S so that essentially all space 
was spanned. Using the expression for H(z) with the two optimized parame- 
ters, the break frequency for the all-pass filter was optimized In a 
similar fashion. The performance criterion used was the mean absolute 
error between the frequency domain phase curve for H(s) and that of 
H(z). 

The frequency characteristics of the final design are shown in 
Figures 3.5 and 3.6. Figure 3.5 shows the magnitudes of H(s) and H(z) 
plotted versus the real frequency w. The maximum magnitude error is 
approximately 1.28 dB which occurs at 20 rad/?.ec., and the mean absolute 
magnitude error is .342 dB. Figure 3.6 shows the phase shift curve for 
H(s) and the phase shift curve for H(z) with and without the all-pass 
filter. For the curve with the all-pass filter over the frequency range 
0-16 rad/sec. the following error information was obtained. The maximum 
phase error between H(s) and H(z) was 12.74 degrees, occurring at 16 
rad/sec., and the mean absolute phase error was 3.86 degrees. Over 
the full frequency range 0-20 rad/sec., the maximum error was 42.4 
degrees, occurring at 20 rad/sec., and the mean absolute phase error 
was 8.61 degrees. 

The optimization procedure was carried out again; this time, using 
the mean squared error criterion on the magnitude curve and the mean 
absolute error on the phase curve. The maximum magnitude error was 
1.17 dB, and the mean absolute macnitude error was .364 dB. Over the 
frequency range 0-16 rad/sec., the maximum and mean absolute phase errors 
were 12.34 degrees and 3.52 degrees, respectively. For the frequency 
range 0-20 rad/sec. the maximum and mean absolute phase error were 42.4 
degrees and 8.30 degrees, respectively. In each case the maximum error 
occurred at the highest frequency in the range specified. 

Trying to use the mean squared error criterion to optimize the 
location of the all-pass filter produced poor results when applied to 
this system. With an all-pass filter, this design procedure produces a 
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phase shift of -ISO degrees at u /2, regardless of the phase shift of H(s). 

s 

For the continuous autopilot the phase shift at 20 rad/sec. is >137.6 
degrees. With the sampling frequency fixed, nothing can be done to reduce 
this large error at ui /2. Since the mean squared error criterion tends 
to accentuate large errors, applying this criterion to the phase curve 
results in increased errors at low frequencies without substantially 
reducing the errors near Wg/2. 

For the optimization schemes tried it appears that for this system the 
mean squared magnitude error and mean absolute phase error criteria pro- 
duced slightly smaller errors than the mean absolute magnitude and phase 
error criteria. The transfer functions 'or the two realizations are 
shown below. In each cas3 they are of the form: 


H(z) 


KN(z) 

nsTiy 


N(Z) = Hj(z) • H 2 (z) • H 3 (z) • H 4 (z) • HgCz) 


D(z) = Hs(z) • Hg(z) • Hyfz) • He(z) 


In the table below the heading MAM/M'\P indicates the mean absolute 
magnitude and mean absolute phase error criteria, and MSM/MAP indicates 
the mean squared magnitude and mean absolute phase error criteria. 



MAM/MAP 

MSM/MAP 

K 

7.80489 E-02 

7.73779 E-02 

Hi(z) 

z - .769409 

z - .769409 

H2(z) 

z2 - I.639I8Z + .69577 

z2 - I.639I8Z + .69577 

H3(z) 

z2 - .223826Z + .4308 

z2 - .223826Z + .4308 

H4(z) 

(z + .230133)2 

(z + .223527)2 

Hs(z) 

(z - .907063)2 

(z - .907063)2 

H6<z) 

z - .470564 

z - .486006 

H7(z) 

z2 - .750535Z + .276885 

z2 - .750535Z + .276885 

He(z) 

z2 - .280832? + .191517 

z2 - .280332Z + .191517 

HgCz) 

(z + 2.72379)/(z + .367135) 

(z + 2.12yi9)f{z + .367135) 


51 



H 4 (z) represents the second order real zero added during design, Hg(z) 
is the single real pole of H(s) whose location in H(z) was optimized; 
and H 9 (z) Is the all-'pass filter. 


Design Considerations 

Several factors concerning this design procedure should be mentioned. 
First, when H(z) is expressed as a ratio of polynomials, the degree of 
the numerator and denominator in powers of z will be equal. The degree 
of the polynomials In z will be equal to one, plus the highest degree of 
the polynomials in H(s). The additional power of z is due to the all- 
pass filter. In this example H(s) was of seventh degree in the denominator 
and fifth degree in the numerator. H(z) was of eighth degree in both 
numerator and denominator. Because of the equal degrees of the numerator 
and denominator In H(z), the phase shift at w /2 will be 0 degrees if no 
all-pass filter is used, and -180 degrees if a phase lag all-pass filter 
is included. A phase lead all-pass filter, (A + p)/(A - p), produces 
an unstable condition. In the frequency domain a good fit can be obtained 
for the magnitude curve over the range 0 < » < ia^/2. For the phase 
curve the closeness of fit depends on the phase characteristics of the 
continuous system. 


Time Domain Evaluation 

In the following discussion concerning time domain performance the 
H(z) obtained from the mean absolute magnitude and phase error criteria 
will be used. This H(z) will be referred to as the DISCRETE approximation 
to H(s). 

To obtain the time domain response of the continuous time autopilot 
a fourth order Runge-Kutta-Gi 1 1 numerical integration was performed. The 
integration step size was chosen to be T/IO. This step size provides 
approximately 48 samples per period of the highest damped natural fre- 
quency in H(s) and 12 samples per time constant for the shortest expo- 
nential time constant. The inputs chosen were the step function and 
sine waves of each of the integer radian frequencies from I to 20 rad/sec., 
i nclusive. 
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The time response of the Tustin approximation to the autopilot was 
obtained by expressing H(s) as a ratio of polynomials and making the 
following substitution: 



r T (z + I) 
2 iz - \) 




K 


The resulting expression was put In the form of a recursive equation 
and solved for the same set of Inputs previously mentioned. 

The H(z) obtained from the design procedure described was also 
expressed In recursive form and solved for the time response. For 
each of the three sets of solutions the outputs were calculated over 
the time range 0-100 T. For the Runge-Kutta solutions this means that 
ten Integrations will be performed per output sample. 

Using the Runge-Kutta solution as a reference, the mean squared 
errors of the Tustin and Discrete approximations were calculated for 
each of the inputs and averaged over the sinusoidal inputs. Figure 3.7 
shows the errors as a function of input sinusoidal frequency. Each of 
the approximations had negligible steady state error for the step 
input. The table below lists the error data for the two approximations. 
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MEAN SQUARE ERROR 


RAD/SEC 

TUSTIN 

DISCRETE 

1 

.16553969 E-05 

.24048390 E-03 

2 

.28388786 E-05 

.16118844 E-03 

3 

.22027003 E-04 

.16690678 E-02 

4 

.22509128 E-03 

.2850k575 E-02 

5 

.84790011 E-03 

.29299688 E-02 

6 

.15901896 E-Oi 

.17630448 E-02 

7 

.17861110 E-02 

.57113966 E-03 

8 

.14286007 E-02 

.78087308 E-04 

9 

.10556364 E-02 

.46068448 E-04 

10 

.10220778 E-02 

.94572395 E-04 

II 

.13234237 E-02 

.95276045 E-04 

12 

.17965967 E-02 

.74684822 E-04 

13 

.22863912 £-02 

.76157714 E-04 

14 

.26955963 E-02 

.11930201 E-03 

15 

.29781393 E-02 

.20891043 E-03 

16 

.31201926 E-02 

.34898954 E-03 

17 

.31259800 E-02 

.54962620 E-03 

18 

.30092672 E-02 

.82658841 E-03 

19 

.27900444 E-02 

.11962247 E-02 

2C 

.22923737 E-02 

.22923735 E-02 


Averaging over the 20 sinusoidal inputs yields the following data: 


Tustin Error 
Discrete Error 


.16700066 E-02 
.83237062 E-03 


For the step input the errors are: 


Tustin Error 
D'screte Error 


.49399356 E-03 
.99913922 E-03 
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For the sinusoidal Inputs the average mean squared error for the 
Discrete approximation is one half that of Tustin's approximation. The 
greatest inprovements over Tust I n*s. method lie In the frequency range 
8>I7 rao/sec. In thl« range the magnitude curve of the DISO^ETE syston 
fits that of the cctlnuoiis syston very closely, with a maximum difference 
of .34 dB. The phase curves also have a good fit over most of this 
frequency range, with less than 10-degree error from 8 to 15 rad/sec. 
and 18-degree error at 17 rad/sec. 

Over this Sivne frequency range, the Tustin frequency domain magnitude 
error varies from 2.1 to i5.7 dB, referenced to the exact magnitude func- 
tion. The phase shift error in this frequency range averages 16 degrees, 
with a maximum of 32 degrees at i7 rad/sec. and is 25 degrees at 15 rad/sec. 

Conclusion 

A procedure has been presented by which a discrete time system can 
be designed to have frequency domain characteristics simiiar to that of 
a continuous time syst«n. For the magnitude curve a close fit can be 
obtained over the frequency range 0 to »^/2. For the phase curve the 
closeness of fit depends on the phase characteristics of the continuous 
system. The design procedure can be carried out either graphically 
or by computer. The location of all critical frequencies In H(p) and, 
thus, the form of H(z) could be optimized, hssed on a number of different 
performance criteria. An in^rovefsent in time domain performance in the 
middle and upper frequcricies was obtained, compared with Tustin's method, 
with a degredaTion of performance in the low frequencies. 

It appears that frequency don»in methods, and the particular procedure 
described here, are valid design approaches. The characteristics of the 
continuous time systsn being modeled and the input frequency ranges of 
interest will determine which approach is best. 
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IV. SUBSTITUTIONAL METHODS 


An investigation Into the suitability of various substitutional 
formulas and techniques for real time simulation is being conducted. 

The study is divided into two phases* the first being concerned with 
linear systems and the second with nonlinear systems. Results from the 
first phase were reported in the S^i-Annual Report of September I97S 
and will not be repeated herein. Progress to date in the second phase 
of this work is the subject of the remainder of this section. 

Several methods for digital simulation of linear systems have been 
investigated* and results show that the IBM method is probably the best 
in terms of error and computation time. Present efforts are being 
directed at the study of nonlinear systems. 

Along with the substitutional methods* other methods* such as IBM* 
Optimum Discrete Approximation and Discrete Compensation* are being 
considered. There are no other methods currently available from litera- 
ture and publications. 

Comprehensive studies of simulations of nonlinear systems* such as 
error analysis* selection of simulation increment* etc. is currently 
beyond the state-of-the-art. This is because nonlinear syst^ns are 
difficult to classify. Comparisons of these digital simulation methods* 
therefore* are largely experimental. Several systems with different 
degrees of complexity and nonlinearities will be studied before any 
conclusions are reached. 

A simple and ’’slightly nonlinear” system was Investigated* and the 
results showed little effect by the nonlinearity on the overall system 
response obtained with different methods. Another system Is currently 
under investigation and will be simulated as soon as the design is 
comp I eted . 

The following paragraphs suggest the methodology being used to 
develop the necessary computer software for analyzing various procedures. 

I. IBM method: The objective is to derive transfer functions of 

Fig. 4.1 such that the poles of the individual transfer functions (Gjiz) 
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and 62(2)), as well as the poles of the overall closed loop transfer func- 
tion, are correct. This is done by matching the closed-loop eigenvalues 
and the static gains of the continuous and discrete system models. 

Poles of G|(z) and 62(2) are in the same location as GjCs) and 
G2 (s) In the s-plane; therefore, the transient response will be correct. 
This characteristic is not adequately preserved in other substitutional 
methods (Tustin, Boxer-Thaler, Madwed, etc.). Note that G|(z) and G2<z) 
are not the Z-transforms of Gj(s) and G2(s) but are transfer functions 
to simulate the continuous response. 

The design procedure is as follows: 

(a) Replace the continuous transfer functions by the discrete 
transfer function, using the Z-transform. A single-period 
delay is inserted in the feedback paths to insure realizability. 
This delay will be compensated for in the final model. 

(b) For each transfer function the static gains are matched 
between the discrete and continuous blocks. The final 
valu3 theorem is applied to find the gain necessary to 
equate the static gains between the two transfer functions. 

(c) Each nonlinear element is replaced by a nominal gain, and 
the closed loop eigenvalues of the continuous and the 
discrete syst«ns are made equal by insert ng a gain in the 
forward loop. This is the most tedious task for a complex 
system of order higher than 4 and with a multiple input /output. 
Fortunately, this mechanism can be implemented on a digital 
computer, as proposed in an IBM report ( Numerical Techniques 
for Real-Time Digital Flight Simulation ). This program will 

be discussed in more detail In a later part of this section. 

(d) Finally, the steady-state gains of the over-all closed loop 
system are matched. The discrete system is also matched to 
either a specific input or Its approximation. An input 
transfer function will be attached in front of the discrete 
system. There are two ways of approximating the input: 
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one by a stair-step function (zero order hold); and one 
by straight line segments (first order hold). 


Input approximation by a stair-step function ; To find the 
Input transfer function, equate the Z-transform of the over- 
all continuous system with the transfer function of the 
over-all discrete system. Since the zero-order hold will 
introduce a half-period lag, the stair-step approximation 
must be designed with a half-period advance so as to obtain 
the correct simulation. 

Input approximation by straight-line segments : The same pro- 

cedure is applied as above, but no conpensation for the time 
shift is necessary, A sumnary in block diagran form is shown 
in Fig. 4.2. 

Programming details for matching eigenvalues in the IBM method are 
as follows. As mentioned earlier, matching eigenvalues of complex con- 
tinuous and discrete systems is a formidable task which cannot be per- 
formed by hand. It has been proposed that computer programs be used to 
plot root loci of the systems and, thereby, determine the gain required 
to match their pole locations. 

There are two possible approaches for using a digital computer to 
compute root loci. One is to determine the characteristic equation as 
a function of loop gain and solve it for various values of loop gain. 

In a multiple loop complex system this is not always possible. There- 
fore, It is desirable to have the computer derive the characteristic 
equation. Another approac'" is to apply Evan’s rule directly with the 
aid of the computer. Thi, method, however, does not fully take advantage 
of the high speed computer. 

IBM has developed a method based on Evan’s rule and used the com- 
puter efficiently at the same time. The program developed by IBM allows 
the user to find root loci of systems (discrete or continuous) directly 
from the block diagram. In the very near future a similar program which 
can be used on our computer will be developed so that more complex systems 
can be simulated. 
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Figure 4.2 Block Diagram of Design Method 
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The principle of the program is as follows. The characteristic 
equation of the system shown in Fig. 4.3 is: 

KG(s)H(s) +1=0 (4.1) 

or 

G(s)H(s) = -l/K (4.2) 

Equation (4.1) is satisfied only if 6(s)H(s) is real, since K is real. 
Furthermore, if a complex number S = o + jw is substituted into G(s)H(s), 
the result will be another complex number q = u + iv. Therefore, the 
solution of Eq. (4.1) is all s that will give v = 0 and K = -1/U. It 
can also be shown that for any path crossing a locus, u wi 1 1 change sign. 
Therefore, the left-half plane of s is scanned until a change In sign of 
V is observed, then the point lying on the locus can be found by iteration. 
The same principle can be applied for a discrete system to find Its 
root locus. 

It has been reported that the developers of the IBM method success- 
fully simulated in real-time complex, six-degree-of-f reedom space vehicles 
and aircraft with considerable improvements over other methods. For 
example, for the same degree of accuracy, this method reduces the compu- 
tation time by 90 to 95$ over the Runge-Kutta method In several cases. 

2. Optimum Discrete Approximation: 

The objective of this design technique is to approximate each elanent 
of the transfer function so that the summation of squared error is minimized. 
This technique eliminates the need for inserting a period delay in the feed- 
back path, since the realizability has been taken into account during the 
design process of the discrete transfer functio-. . Each discrete transfer 
function is considered to consist of a tanuem connection of two elements, 
Fl(z) and F 2 (z). F 2 (z) is called thf» +lxed portion, and it can take the 

form 

F2(z) = Z"P (4.3) 

where p = 0,1. Whether or not F-fz) = I or F 2 (z) = Z"^ depends on the 
system under consideration. Transfer functions where the degree of the 
numerator is lower than the degree of the denominator are called closed- 
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loop realizable to ensure the realizability of the whole system. There- 
fore, in designing the discrete transfer functions we initially select 
p 2 (z) =1. If none of the discrete transfer functions turn out to be 
closed-loop realizable, then we pick the simplest one to be redesigned 
with p 2 (z) = Z“*. This will greatly simplify the work involved. 

Sage and Burt have shown that the discrete ''optimum'* transfer func- 
tion is 


Plopt(Z) 


p2(Z“hU(Z-MA(Z) ^ 
CP2(Z)P2(Z'1)U(Z)U(Z‘^0_ J 
CP2(Z)P2(Z‘MU(Z)U(Z“M]^ 


(4.4) 


with 

A(z) = Z[6 (s)U(s)3 (4.5) 

where U(t) is a test signal and where denotes all poles and zeros 
within the unit circle, denotes all poles and zeros outside the 
unit circle, and 'S the portion with poles within the unit circle. 

The procedure can be surmnarized as follows: 

1. Select a test signal, U(s), usually either a unit ramp or a 
unit step. 

2. Design each discrete transfer function, using the F^opt(z) 
equation. Let P 2 (z) = I in each case; and, if Piopt(z) is not closed- 
loop realizable, then redesign it with p 2 (z) = Z”^. Then 


G(z) = Piopt(z)p 2 (z). (4.6) 

The procedure above does not take into account the nonlinear characteristic 
of the system. It will work adequately for a slight nonlinearity. If 
the system is decidedly nonlinear, the approximation error can be reduced if 
we make use of the fact that the system is actually nonlinear in determi- 
ning the "optimum" discrete approximation. We can expand the method 
above, using gain parameters before and after the nonlinearity. These 
gain parameters will be adjusted during the simulation process. 
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3. Discrete Ckjtnpensation: 

In substitutional methods it is necessary to insert a period delay 
in the feedback path to ensure realizability. This may shift the pole 
locations in some systems. The IBM and Optimum Discrete Approximatior 
methods compensate this delay In the design process. In the Discrete 
Compensation method the closed-loop transfer function is discretized, 
using any discrete integrator, but improvements are made through a 
compensator to adjust the eigenvalues of the closed- loop system. 


For example, take a simple system, such as shown in Fig. 4.1 and 


let 


Gi(z) = Gi(s)i 


= some i ntegrator 


(4.7) 


62(2) = 62(5) 


~ = some I ntegrator 


(4.8) 


H(z) = H(s) 


(4.9) 


= some i ntegrator 


We then discretize the whole system as 

G 3 G 1 ( S )G 2 (s) 

" I + H(s)53Gi(s)G2(s) 


(4.10) 


— = some i ntegrator 


where G 3 is the nominal gain of the nonlinear element. Insert a compensa- 
tor after the input of the discrete system and compare it with G^(z) to 
obtain the coefficients of the compensator, whose form is 


bo + biz + + b z*” 

D(z) = 2— (m < n) (4.11) 

»0 * =1^ * * V 

Gj ( z)G 2 ( 2)63 

G (z) = D(z) ; (4.12) 

I + G3Gi(z)G2(z) 

The final form is shown in Fig. 4.4. 
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Figure 4.4 Final Form of Discrete Compensation Method 




If further accuracy Is required (at the expense of more computation 
time), the coefficients of D(z) can be adjusted, depending on the In- 
stantaneous gain of the nonlinear element; that Is, Instead of using 
G3 as a nominal gain to obtain fixed coefficients of 0(z), we can make 
the coefficients of D(z) a function of §3, where G3 can be determined 
by direct measurement by Interpolation. 

Summary 

In the three methods considered above the Discrete Compensation Method 
is the best In terms of design efforts. Its drawback Is excess computa- 
tion time, especially with the coefficient adjustor. Once the computer 
program for the IBM method is completed, the IBM method may be the logical 
choice, even though the Optimum Discrete Approximation may be more accurate 
for decidedly nonlinear systems. However, no final conclusions can be 
drawn at this time. 

Future Tasks 

• Developing a computer program to aid the design of the IBM method. 

• Further search and studies of Sage's method, using calculus of 
variations and other techniques for a nonlinear system. 

• Simulations of various nonlinear systems with special attention 
given to the criteria mentioned above, observing possible effects 
of a certain nonlinearity on a certain method. 
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V. USE OF FADE APPROX I MANTS TO THE f^ATRIX EXPONENTIAL FOR 
COMPUTER SOLUTIONS OF STATE EQUATIONS 

Introduction 

In system theory, there Is a large class of problems which may be 
phrased In terms of linear time- Invar I ant differential equations, and 
which lend themselves to straightforward solutions. A second class, at 
the opposite end of the spectrum, consists of problems characterized by 
behavior which Includes large time variation and s-*rong non-l Inearltles. 
There is a middle ground, however, consisting of a class of problems In 
which the time-dependent parameters vary relatively slowly, with weak 
non-l inearltles. In these, it may be necessary (or merely desirable) 
to Include the effects of time- and state-dependent variables, and at the 
same time undesirable to utilize algorithms used normally on systems of 
large computational complexity. The results of the theory of linear time- 
invariant systems may be used, with proper modification, to approximate 


this class of systems very closely. 

The matrix differential equation 

i = Ax + Bu (5.1) 

where A and B are constants, has the well known solution 

x(t + T) = G(T)x(t) + H(T)u(t), t = kT, k = 0,.,2, ... , (5.2) 

where T is taken to satisfy the Nyquist criterion on u, 

G(T) = exp(AT) = I + AT + (A^T^/z!) + ... , (5.5) 

and 

H(T) = exp(Ax)dTB = [(IT + (AT2/2!) + (A^tVS!) + ...)]B. (5.4) 

The last term may be reduced to 

H(T) = A-*[(exp(AT) -I)]B, (5.5) 


If the inverse of A exists. The above may be derived directly from the 
Taylor series of x(t + T). 
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A more general statement is true where A - A(x,t) and Is a slowly 
var, ..g function of both x and t. For the class of problems where the 
system ti mo-dependence and non-linearities are small, we may still use 
the matrix exponentail to approximate the system*s behavior C5.G. That 
•s, 

x(t + T) = exp(A(x(T),t)T)x(t) + exp(A(x(T),t)T)dTBu(t). (5.6) 

Allowing A to vary with x and t, however, suggests consideration of 
effici<;nt techniques for computation of exp(AT). In some cases, simple 
summation of the terms of the power series representation will be suffi- 
cient since A may vary so slowly that only occasional updatings of exp(AT) 
are necessary. In other cases, it tnay be necessary to recompute exp(AT) 
at each sampie point to achieve the desired accuracy. The Fade approx I mants 
to exp(AT) wnich will be considered offer, in some types of problems, 
significant advantages in computational speed and accuracy over the power 
series representation. 


Definition of Fade Approx i mants 

A Fade approximant to o scalar power series F(z) is a ratio of two 

polynomials P(z) and 0(z), of erder p and q respectively, abbreviated 

(p,q). Its significant feature is that the power series expansion of 

(p,q) is identical with that of F(z) up to and including the coefficient 

N 

of z , where N = p + q is the order of the approximant. There are N + I 
Fade apprexi mants of order N, and they are unique C5.2]. For example, 
the three approx i mants for N = 2 and F(z) = exp(z) are: 

(2, 0) = I + z + zV2 

(I, !)=(!+ z/2)/(l - z/2); lz| < 2 

(0, 2) = 1/(1 - z + zV2) lz| < 1.414 (5.7) 

Similar statements may be made of Fade approxlmants to s matrix power 
series; for example, the (I, I) approximant tc ‘■xp(AT) would be: 

(I, I) = (I - AT/2)" (I + AT/2) = (I + AT/2) (I - AT/2)" (5.8) 
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Of the general class of Fade approximants, some are more effective 
computationally than others. In general it requires no more computation 
to calculate a (P, P) Fade approximant than a (P-M,P) or a (P, P-M), 
where M < P. Because the (P, P) approximant Is accurate through more 
terms, it is the most beneficial form to use. In this paper we shall 
compare (I) the truncated series (N, 0), and, (2) those in which p = q, 
i.e., (I, I), (2, 2), (3, 3), and so forth. The scalar approximants 
(2, 2) and (3, 3) for exp(z) are: 

(2, 2) = (I + z/2 + zVl2)/(l - z/2 + zVl2); \z\ < 3.464 (5.9) 

(3, 3) = (I + z/2 + zViO + zVl20)/(l - z/2 + zViO - zVl20); 

|z| < 4.644 (5.10) 

It is shown in Appendix A that the error involved in using a Fade 
approximant to the matrix series exp(AT) is a linear function of the 
error in approximating the scalar series exp(X^T) where is the 
magnitude of the maximum eigenvalue of A. This allows discussion of 
accuracy in terms of the scalar series exp(X^T), with results which 
carry over directly to the matrix series approximation. 

Accuracy vs. Required Confutation 

The two classes of approximations for the exponential were compared 
directly, using three criteria: (I) the numer of matrix operations re- 

quired for the approximation to exp(AT), (2) the percent error in terms 
of the scalar approximation of exp(X^T), and (3) the allowable range of 

X T for which reasonable results could be achieved. It will be shown 
m 

that the inclusion of higher-ordered terms in the (p, p) type of approximant 
reduces error, extends the range of X^T, and is in general more efficient 
than the corresponding (p + p, 0) truncated series. 

A matrix operation may be defined as a multiplication or an Inverse; 
each requires about n^ multiplications and divisions, where n is the order 
of the matrix C5.3]. This is an effective standard, since the matrix 
manipulations dominate the calculation time. By this standard, the (N, 0) 
series requires N-l matrix operations, while the (p, p) series requires 
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p ♦ I operations. For example, the (2, 2) and the (4, 0) both require 
three matrix operations; however the (3, 3) uses only four while the 
(6, 0) requires five operations. For systems of large order, this may 
be significant. 

A computer was used to calculate % error in the various forms of 
approximations to exp Some of the more Important features are 

demonstrated in Figs. 5.1 and 5.2. For example, in Fig. 5.1, for X^T = .p, 
the (4, 0) truncation yields about .04? error. The (2, 2) gives approxi- 
mately .004? error, or about one order of magnitude improvement. For 

large values of X T the results are even more dramatic. For X T = 2.0, 
m m 

a (9, 0) approximation gives .17? error con^)ared to the .15? error of the 
(3, 3) approximation, which also affords a savings of four matrix opera- 
tions. This alone may result in a considerable reduction in computer 
time. Fig. 5.2 shows the maximum allowable values of |x T| for ^iven 
amounts of error. The (p, p) type of approximant is seen to have consistent- 
ly greater range than its corresponding truncated series. 

Previously, it has been necessary to choose to be relatively 

small in order to achieve the required accuracy with a short series 
approximation. Using the Fade (p, p) approximants, there is much more 
room for choice. Larger system eigenvalues may be included in the 
problem statement, or a longer sample time T chosen (subject to the 
Nyquist rate on u(t)), and the desired accuracy may still be achieved, 
or bettered, within the context of a real-tiine computer simulation. 

Summary 

In summary then, it has been shown that thei (p, p) Fade approximant 
to the matrix exponential exp (AT) is a useful tool in control system 
simulation and operation. It has the advantages of greater accuracy, 
larger range, more flexibility, and in s«ne cases greater computational 
efficiency than the truncated series approximation. 
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Appendix A 

Relation' Between ^tet^ix Series and Scalar Series 

To obtain a measure of the error introduced by using approximations 
for the matrix exponential, one can transform to the normal coordinates, 

i .e. 

X = s“^;s 
and 

X = S X . 

It will be assianed that A has distinct eigenvalues. This is not 
overly restrictive since any matrix with repeated eigenvalues may be 
approximated arbitrarily closely by one with distinct eigenvalues C3.4D. 

The matrix exponential operating on a vector may be written as 
exp(AT)x(t) = SS‘^exp(AT)SS’lx(t) 

= SS**(I + AT + i A^t 2 + )SS”^x(t) 

= S{S"MS + S"*AST + j S"^ASS-*AST2 + )S“lx(t) 

= S(l + XT +|-X2 t 2 + )i(t) 

Likewise one can write the (I, I) Fade approximants as 
(I -^AT)”^(I + jAT)x(t) = SS“^I -iAT)"l(l + j AT)SS"^x(t) 

= SS*^I +iAT + jA2T2 + )SS"Ml + i AT)SS-^x(t) 

= S(S-^IS + j S"^AST + I S‘^SS~^AST2 + HS-^IS + j S"^AST)S'*x(t' 

= S(l + 7 XT + j X^t2 + )(| + j XT)^(t) 
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The other ordered Fade approximants follow In a like manner. In 
each case the approximant as well as the exact expression for exp(AT) is 
a diagonal matrix in the normal coordinates. Also each diagonal element 
is a function of one and only one eigenvalue. This fact makes it possible 
to compare different approximations to exp (AT) in the normal coordinates on 
an element by element basis. Also, although we have not proved it mathe- 
matically, our numerical results indicate that the maximum error occurs in 
that element corresponding to the eigenvalue of maximum value. Thus 
one can test for the accuracy of matrix polynomial representations of 
exp(AT) by testing the corresponding scalar polynomials and utilizing 
the eigenvalue of the A matrix having greatest magnitude. 
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Appendix B 

Fade Apptioximants for H(T) 

One further use for the Fade approximation may be considered here. 
It was shown above that If A“^ exists, then 

H(T) = A"iCexp(AT) - I3B. 

The A matrix may be singular, however, and no convenient closed-form 
expression can be used for H(T). For example, a zero eigenvalue in A 
will impose this restriction. However, H(T) may be written 

H(T) = T( I + AT/2 + A^tVS! + )B. 

This series written in scalar form is 

I + w/2 + w^/3! + — 

where w = X T. This series has the scalar Fade approximants: 
m • 

(I, I) = (I + w/6)/(l - w/3); |w] < 3 

(2, 2) = (I + w/IO + wV60)/(l - 2w/5 + wV20); |w| < 4.472 

(3, 3) = (I + w/14 + wV42 + wV840)/(l - 3w/7 + w^/14 - wV2IO); 

|w| < 5.649 

Use of these approximants parallels the discussion above for the 
matr i x exponent i a 1 . 
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VI. ANALYTICAL INTEGRATION OF STATE EQUATIONS, 

USING AN INTERPOLATED INPUT 

One approach to solving differential equations discretely has been 
to use the state equations, along with polynornial appro I iiat Ions to the 
input signal. The coefficients of the polynomial are determined by the 
value of +he Input signal at the sample times. Once the oolynomial 
is chosen, the state equations can be integrated analytically. One 
important feature of this approach is that the system is modeled exactly. 
The approximation is in sampling the input. The various polynomial fits 
to the Input serve as Interpolators. 

Using a first order polynomial to represent the input, one finds 

that 

V - ^ 4 . fNT A(NT-t)Q,„., .. 

^ ~ ® ^-1 * ^(N-I)T ® BU(t)dt (6.1) 


can be approximated by 

( 6 . 2 ) 

Expressions for the zeroth order polynomial fit and the 2nd order 
polynomial fit have been determined and are reported along with more 
details on the method in our Semi-Annual Report of September, 1975, on 
this same project. University Report No. EE-4041-1 01 -75. 

Because our efforts have been directed in other areas, we have not 
yet obtained test results on the accuracy of this technique. 
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